Capillary Rise in Nanopores: Molecular Dynamics 
Evidence for the Lucas- Washburn Equation 
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When a capillary is inserted into a liquid, the liquid will rapidly flow into it. This phenomenon, 
well studied and understood on the macroscale, is investigated by Molecular Dynamics simulations 
for coarse-grained models of nanotubes. Both a simple Lennard- Jones fluid and a model for a 
polymer melt are considered. In both cases after a transient period (of a few nanoseconds) the 
meniscus rises according to a V time- law. For the polymer melt, however, we find that the capillary 
flow exhibits a slip length S, comparable in size with the nanotube radius R. We show that a 
consistent description of the imbibition process in nanotubes is only possible upon modification of 
the Lucas- Washburn law which takes explicitly into account the slip length 8. We also demonstrate 
that the velocity field of the rising fluid close to the interface is not a simple diffusive spreading. 

PACS numbers: 47.11. +j, 47.55-t, 47.60.+i, 66.30.-h 



Introduction Understanding fluid flow on the 
nanoscale [l[ is crucial for modern developments of 
nanotechnology like the "lab on a chip" and related 
nanofluidic devices, as well as for various applications of 
porous materials, fluid flow through pores in biomem- 
branes, etc. A key process is the ability of fluids to 
penetrate into fine pores with wettable walls. Filling hol- 
low carbon nanotubes or alumina nanopore arrays with 
chosen materials opens exciting possibilities to generate 
nearly one-dimensional nanostructures 0> 0]. In this 
context, also the filling of silicon dioxide nanochanncls 
[3] and of rolled- up InAs/GaAs tubes [f| has found 
great interest. Related fluid filling phenomena occur 
when viscous fluid fronts propagate into porous media 
by spontaneous imbibition [6]. On macroscopic scales, 
a basic understanding of such capillary rise processes 
exists for almost a century 0, H, @, 13, hjj, Uj- How- 
ever, the applicability of the resulting concepts on the 
nanoscale has been the subject of a recent controversy 
E3, EH . In particular, the conditions under which 
the Lucas- Washburn equation 0, @] holds are debated. 
This equation predicts a \ft law for the rise of the fluid 
meniscus H (t) in the capillary with time t, 



H(t)=^ LV ^ 0Sd ) 1/2 Vi. 



(1) 



Here jlv is the surface tension of the liquid, r] its shear 
viscosity, R the pore radius, and 9 the contact angle be- 
tween the meniscus and the wall. Eq.([T]) follows by in- 
tegration of the differential equation, describing steady 
state flow, where the capillary force 2jlv cos(9)/R is bal- 
anced by the viscous drag ind(H (t) / R) 2 / dt and one as- 
sumes that any possible slip length 8 <C R- Of course, 



Eq. (TTJ) cannot be true for t — *■ 0, but can hold only af- 
L, = 80 




FIG. 1: Snapshot of fluid imbibition in a capillary at time 
t = 1300 MD time steps after the onset of the process. Fluid 
atoms are shown in blue, those of the precursor film - in light 
blue. The tube wall is shown in red, and the atoms of the 
reservoir adhesive wall - in yellow. For further explanations 
see text. 

ter a (nanoscopically small) transient time (Zhmud et al. 
[lH suggest an initial behavior H(t) cx t 2 when the liquid 
is accelerated by the capillary forces). However, Mastic 
ct al. [IH find H(t) rising slower than linear with time, 
even for t sa Ins, from Molecular Dynamics (MD) sim- 
ulation of a simple Lennard-Jones (LJ) fluid. They sug- 
gest to slightly correct Eq. ([1]) , replacing 9 by a dynamic 
contact angle 9{t). In contrast, a study of a model for 
decane molecules in a carbon nanotube [3, EH yielded 
a simple linear behavior H (t) cx t over a wide range of 
times, leading to the conclusion that filling of nanotubes 
by fluids does not obey the Lucas- Washburn equation. 
Experiments so far are inconclusive on this issue, since 
the existing work [3, E3] deals only with pores that are 
at least 1 fim wide. Moreover, in narrow nanotubes an 
eventual slip at the hydrodynamic boundaries might af- 
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feet the balance of forces by reducing the viscous drag at 
the tube wall. 

The aim of the present letter is to help clarifying the 
problem of capillary filling in narrow nanotubes. We 
present simulations of a generic model, varying both the 
fluid-wall interaction and the nature of the fluid (sim- 
ple LJ particles vs. melt of short polymer chains, re- 
spectively). Providing independent estimates for all the 
parameters entering Eq. |T]), we are able to perform a 
decisive test of Eq. fTJ). Since a fluid flowing into a 
capillary is a nonequilibrium process, we avoid use of 
a strictly "microcanonical protocol" of our MD simula- 
tions, unlike [3, 15 1. Using a dissipative particle dy- 

which does not disturb 



namics (DPD) thermostat [if 



the hydrodynamic interactions due to its Galilean invari- 
ance, we maintain strict isothermal conditions, in spite of 
the heat produced due to the friction of the flowing fluid. 
In the real system, the walls of the nanofluidic device 
would achieve the thermostating, of course. 

Model description - The snapshot picture, Fig. Q] il- 
lustrates our simulation geometry We consider a cylin- 
drical nanotube of radius R = 10, whereby the cap- 
illary walls are represented by atoms forming a trian- 
gular lattice with lattice constant 1.0 in units of the 
liquid atom diameter a. The wall atoms may fluctu- 
ate around their equilibrium positions at R + a, sub- 
ject to a finitely extensible non-linear elastic (FENE) po- 
tential Ufene = -15e w J?o m (l — r2 /Ro) > = 1.5. 
Here e w = l.OfcsT, ks denotes the Boltzmann con- 
stant, and T is the temperature of the system. In 
addition, the wall atoms interact by a LJ potential, 
ULj(r) = 4e tuu) [(a ww /r) 12 - (a ww /r) 6 ], where e ww = 
1.0 and a ww — 0.8. This choice of interactions guarantees 
no penetration of liquid particles through the wall while 
in the same time the wall atoms mobility corresponds to 
the system temperature. In all our studies we use a cap- 
illary length H rnax = 55. The right end of the capillary is 
closed by a hypothetic impenetrable wall which prevents 
liquid atoms escaping from the tube. At its left end the 
capillary is attached to a rectangular 40 x 40 reservoir for 
the liquid with periodic boundaries perpendicular to the 
tube axis. Although the liquid particles may move freely 
between the reservoir and the capillary tube, initially, 
with the capillary walls being taken distinctly lyophobic, 
these particles stay in the reservoir as a thick liquid film 
which sticks to the reservoir lyophilic right wall. The film 
is in equilibrium with its vapor both in the tube as well 
as in the left part of the reservoir. At a time t = 0, set 
to be the onset of capillary filling, we switch the lyopho- 
bic wall-liquid interactions into lyophilic ones and the 
fluid enters the tube. Then we perform measurements 
of the structural and kinetic properties of the imbibi- 
tion process at equal intervals of time. As a simulation 
algorithm we use the velocity- Verlet algorithm [l9| and 
DPD thermostat [H],[2(| with friction parameter £ = 0.5, 
Heavyside-type weight functions, and a thermostat cutoff 
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FIG. 2: Profiles of the average fluid density p(z) in the cap- 
illary at various times for the case e w i — 1.4, e = 1.4. The 
small oscillations reflect the corrugated structure of the wall 
of the capillary. 



r c = 2.5a. The integration time step St = O.Olto where 
to is our basic time unit, to = \J 'ma 2 / '48fc_BT = 1/V48, 
choosing the particle mass m = 1 and fceT = 1. 

The capillary filling is studied for two basic cases: (i) 
a simple fluid interacting via LJ potential with e = 1.4 
and a = 1.0, and (ii) a non-Newtonian fluid (a polymer 
melt) consisting of short chains of length N — 10. The 
non-bonded interaction is given by a LJ potential with 
e = 1.4 and a = 1.0 whereas the bonded forces between 
chain monomers result from a combination of FENE and 
LJ potentials with e = 1.0 [211 ] . In both cases the liquid- 
wall interaction is given by a LJ potential with strength 
e w i which is varied over a broad range so as to change the 
wetting characteristics of the system. All interactions 
are cut off at r cut — 2.5a. By varying the interaction 
strengths and the thermostat parameters, one can change 
the dynamic properties of the test fluids in a wide range. 
The total number of liquid particles is 25000 while those 
forming the tube are 3243. 

Simulation results - Typical data for the time evolu- 
tion of the advancing front of the LJ fluid penetrating 
into the pore are shown in Fig. [2j Choosing a constant 
time interval (At = 300) between subsequent profiles in 
Fig. El it is already obvious that the interface position 
advances into the capillary slower than linear with time! 
The profiles p(z) at late times become distinctly nonzero 
far ahead of the interface position (near the right wall at 
z = 55 where the capillary ends), due to a fluid mono- 
layer attached to the wall of the capillary: this precursor 
advances faster then the fluid meniscus in the pore cen- 
ter, but also with a y/t law (see below). 

Fig. [3k. shows that the time evolution of the meniscus 
height H (t) depends very sensitively on the strength of 
the wall-fluid interaction (or the contact angle 9, respec- 
tively): for ews — 0.6 and 0.8 only a small number of 
fluid particles enter the capillary (since 8 > 90 for these 
choices (22[). For ews > 1-2, however, there is only a 
short transient up to about t = 250to (i.e., a time still 
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FIG. 3: Position of the liquid meniscus H(t) for a LJ fluid 
(a), and a melt of decanes (b), for various choices of e w i. The 
straight lines in (a) - full line, and in (b) - dashed line, indicate 
the asymptotic law, Eq. ((TJ, in the case of complete wetting. 
The topmost dash-dotted curve indicates H(t) for the pre- 
cursor foot (calculated from the total number of particles in 
the first layer adjacent to the wall of the tube). The initial 
rise of the precursor (for t < lOOto) proceeds much faster, 
however. The inset in (b) shows the radial velocity variation 
in a steady flow regime for two strengths of applied external 
force g — 0.02 and g = 0.05 (see text) clearly indicating a slip 
length 8 « 2.7. 



in the nanosecond range), and then a behavior compat- 
ible with Eq. [T] is verified. The initial deviation from 
this law seen in Fig. [3] is not related to the "dynamic 
contact angle" this would produce a curvature of 

opposite sign in Fig. [3J However, in the marginal case 
£wl = 1.0 a pronounced deviation from Eq. ([T]) occurs: 
this curve could be (approximately) fitted to a linear re- 
lation, H{t) oc t. It is possible that the results of Supple 
and Quirke 



1< 



15| just correspond to such a marginal 
case. Finally we emphasize that even the height of the 
precursor foot (topmost curve in Fig. [3^) advances with 
a H(t) oc \ft law. 

To test whether the capillary rise behavior of polymers 
differs from that of Newtonian fluids, the penetration of 
a melt of short flexible polymers (above model (ii)) is 
shown in Fig.[3jD. But apart from a general slowing down 
(attributable to the higher viscosity of the polymer melt), 
the behavior exhibits the same y/i-lscw. 
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FIG. 4: Velocity field around the moving meniscus for e w i = 
1.4. Velocities are averaged within a two-dimensional grid, 
always fixed at z — to the actual moving meniscus position 
and renormalized according to the current meniscus speed. 
The inset shows the LJ-fluid density profile in the vicinity of 
the meniscus. The interface position is denoted by a yellow 
line. 



We have also tested whether the results are affected by 
possible simulation artefacts due to insufficient thermo- 
stating conditions. In fact, a slower capillary rise occurs 
if one uses a "Langevin thermostat" (i.e., an ordinary 
friction and random noise term acting on all particles), 
which violates hydrodynamics (23|. A similar result ap- 
plies if the wall atoms are rigid rather than mobile [221 ] . 
But even in these cases the data still follows the \/i-l&w, 
and also changing details such as the above parameters 
chosen for the FENE potential of the wall atoms does 
not matter. From the velocity profiles near the moving 
meniscus Fig. ((4} it is evident that care is needed for the 
temperature equilibration of such non-equilibrium MD 
simulations of transient phenomena. The flow is laminar 
behind the interface, parallel to the walls of the capillary, 
with the velocity largest in the tube center and going to 
zero close to the walls. Thus our simple fluid exhibits 
evidently stick boundary conditions. However, in the 
interface the velocity field bends over into a direction 
along the interface, and occasionally particles evaporate 
into the gas region. This flow pattern shows that the 
H(t) oc \ft must not be confused with a simple diffusive 
spreading, of course! 

Comparison with the Lucas- Washburn equation - For 
a test of Eq.[TJ it is crucial to also estimate the prefactor, 
of course, to prove that the \ft growth is not just a mere 
coincidence. For the LJ fluid (at density pi — 0.774) we 
find 77 « 6.34±0.15, for the polymer melt (at pi = 1.043) 
the result is 77 « 205 ± 25. We derived compatible values 
for the viscosity of both fluids also within an equilib- 
rium Molecular Dynamics simulation by using the cor- 
relation function of off-diagonal pressure tensor compo- 
nents and the standard Kubo relation [191. From the 
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flat gas-liquid interface observed in the left part of our 
simulation box (Fig. [T|) we can estimate the surface ten- 
sion j£ V from the anisotropy of the pressure tensor [25[ , 
lew = J dz{p zz (z) - \p xx (z) +p yy (z)}/2}. This yields 
7^ = 0.735 ±0.015 (LJ) and 1.715±0.025 (polymer), re- 
spectively. 

A consistency check of our results with Eq.([T]) is per- 
formed for the case of complete wetting, cos((9) = 1, 
which corresponds to e w i > 1.4 where the data prac- 
tically collapse on a single curve. For the simple LJ- 
fluid with the tube radius R — 10 one obtains a slope 
H(t)/\/t = 0.76 ± 0.02 which agrees perfectly with the 
measured meniscus velocity, cf. Fig. [3^. For the polymer 
melt, in contrast, Eq.([T|) predicts a slope of 0.20 which is 
considerably less than the observed slopes in Fig. [3b. To 
clarify this discrepancy we performed MD simulations 
of steady state Poiseuille flow of identical melt subject 
to external force g comparable to the capillary driving 
force. The radial variation of axial velocity indicates a 
clear slip-flow behavior, cf. inset in Fig. [3Jd, with a slip 
length of S « 2.7 which cannot be neglected when com- 
pared to the tube radius R — 10. The importance of 
slip-length in processes in the nanoscale range has been 
emphasized earlier by Barrat and Bocquet [261 ]. In the 
present case the existence of a slip length S can be eas- 
ily accounted for in the Lucas- Washburn result, Eq.JT]), if 
one notes that, according to the definition of a slip length, 
the drag force under slip-flow conditions in a tube of ra- 
dius R and slip length S is equal to the viscous drag force 
for a no-slip flow in a tube of effective radius R + 5, that 
is, to 4r]d(H (t) / (R + <5)) 2 / dt. In both cases the capillary 
driving force remains unchanged, cos(9) / R. Thus one 
derives a modified Lucas- Washburn relationship: 
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H(t) = 



1lv{R + S) 2 cos9 



-,1/2 



2Rr] 



Vt. 



(2) 



Using Eq.([2]), and the material constants given above, we 
obtain for the slope H(t)/\/t = 0.26 ± 0.02 which agrees 
within errorbars with the observations in Fig. [3b. 

Conclusions In summary, we have shown that ba- 
sic concepts of capillarity such as the Lucas- Washburn 
equation, Eq. [TJ work almost quantitatively even at the 
nanoscale, both for small molecule fluids and complex 
fluids such as short polymer chains. In case of slip- 
flow, however, we suggest a simple modification which 
takes into account the slip length 6. Our new result, Eq. 
©, restores the consistency of the Lucas- Washburn law 
within the framework of the general \ft law even in those 
cases when slip-flow cannot be neglected. 
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